Relativistic Radiation Magnetohydrodynamics in Dynamical Spacetimes: Numerical 

Methods and Tests 
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Many systems of current interest in relativistic astrophysics require a knowledge of radiative 
transfer in a magnetized gas flowing in a strongly-curved, dynamical spacetime. Such systems in- 
clude coalescing compact binaries containing neutron stars or white dwarfs, disks around merging 
black holes, core collapse supernovae, coUapsars, and gamma-ray burst sources. To model these 
phenomena, all of which involve general relativity, radiation (photon and/or neutrino), and magne- 
tohydrodynamics, we have developed a general relativistic code capable of evolving MHD fluids and 
radiation in dynamical spacetimes. Our code solves the coupled Einstein-Maxwell-MHD-Radiation 
system of equations both in axisymmetry and in full 3-1-1 dimensions. We evolve the metric by 
integrating the BSSN equations, and use a conservative, high-resolution shock-capturing scheme to 
evolve both the MHD and radiation moment equations. In this paper, we implement our scheme 
for optically thick gases and grey-body opacities. Our code gives accurate results in a suite of tests 
involving radiating shocks and nonlinear waves propagating in Minkowski spacetime. In addition, 
to test our code's ability to evolve the relativistic radiation-MHD equations in strong-fleld dynam- 
ical spacetimes, we study "thermal Oppenheimer- Snyder collapse" to a black hole, and find good 
agreement between analytic and numerical solutions. 

PACS numbers: 04.25.D-, 04.40.Nr, 47.75. +f, 95.30.Jx 



I. INTRODUCTION 

Many relativistic systems of current astrophysical in- 
terest are characterized by the dynamical coupling of 
strong-field gravitation, high magnetic fields and intense 
radiation (where the latter may be photons or neutrinos). 
Quasars, active galactic nuclei (AGNs), galactic "micro- 
quasars" , core-collapse supernovae, coUapsars, gamma- 
ray burst sources (GRBs), merging neutron star bina- 
ries (NSNSs), merging black hole- neutron star binaries 
(BHNSs) , and merging neutron star- white dwarf binaries 
(NSWDs) are all examples of such systems. Developing 
robust computational methods that can treat simulta- 
neously the different dynamical phenomena that govern 
these systems is necessary in order to simulate their phys- 
ical behavior reliably and identify their observational sig- 
natures. 

Many of the systems listed above involve compact ob- 
jects, such as black holes and neutron stars. Hence gen- 
eral relativity is required to describe their dynamical evo- 
lution accurately. Both observations and theory strongly 
suggest that magnetic fields play an important role in 
many of these systems. For example, magnetic fields are 
crucial in launching jets from black holes in AGNs and 
GRBs (see, e.g., [iIjBIl driving accretion onto black holes 
in disks (see, e.g., 0, |j|), and inducing 'delayed' collapse 
in hypermassive neutron stars that may form following 
NSNS mergers [1, H, 0, Hi- Radiation, apart from its 
role as an observational tracer and diagnostic probe, also 
can play an important dynamical role in many relativistic 
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systems. For example, the role of neutrino transport may 
be essential to understanding core-collapse supernovae 
(see, e.g. @, As a second example, consider that 

the interior pressure of supermassive stars and massive 
Population III stars is dominated by thermal radiation 
pressure. These objects may collapse in the early uni- 
verse to form the seeds of the supermassive black holes 
that reside in the centers of many, and perhaps most, 
galaxies [ll, [3] ■ Radiation thus plays a crucial role in 
determining the onset and dynamics of the collapse of 
these stars and the masses and spins of the black holes 
that are formed As a final example, accre- 

tion onto compact objects leading to outgoing radiation 
near and above the Eddington value is controlled by the 
competition between inward gravitational forces and out- 
ward radiation pressure forces. All of these systems need 
to be handled in a computational scheme designed to 
probe these physical phenomena self-consistently, simul- 
taneously accounting for radiation, magnetic fields and 
relativistic gravitation. 

We have developed previously a robust numerical 
scheme in 3-1-1 dimensions that simultaneously evolves 
the Einstein equations of general relativity for the 
gravitational field (metric), the equations of relativis- 
tic magnetohydrodynamics (MHD) for the matter, and 
Maxwell's equations for a magnetic field Our 
approach is based on the BSSN (Baumgarte-Shapiro- 
Shibata-Nakamura) formalism to treat the gravitational 
field [Hjlill, a high-resolution, shock-capturing (HRSC) 
scheme to handle the fluid and a constrained-transport 
scheme to treat magnetic induction [20]. Our resulting 
GRMHD code has been subjected to a rigorous suite of 
numerical tests to check and calibrate its validity 
We have applied our code to explore a number of dynam- 
ical scenarios, including the collapse of magnetized, dif- 
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ferentially rotating, hypermassive neutron stars to black 
holes 0, [1] , the collapse of rotating stellar cores to neu- 
tron stars the collapse of rotating, supermassive 
stars and massive Pop III stars to black holes [22^, and 
the merger of binary black holes [1^ and binary black 
hole-neutron stars [24]. The purpose of this paper is to 
present a generalization of our current GRMHD scheme 
that accounts for the presence of radiation (photon or 
neutrino). 

Our approach for handling the radiation follows in the 
long tradition of formalisms designed to treat radiation 
transport in the framework of general relativity. How- 
ever, we have developed a new version specifically de- 
signed to fit neatly onto our existing 3-1-1 GRMHD 
scheme. The general relativistic radiative transfer equa- 
tion has been derived in full detail by Lindquist in 
1966 [1^. His treatment has been followed by numerous 
adaptations and implementations in various approxima- 
tions. For example, Thorne has derived a set of radiation 
moment equations to arbitrary order by the technique 
of projected symmetric trace-free (PSTF) tensors 26]. 
So far, most GR radiation hydrodynamics calculations 
(e.g. [23,|2i,[2i,[3S[M[3l[3iJa^ including those based 
on the PSTF scheme (e.g. [Slli, [13, [11, S 113]), have 
been implemented in spherical symmetry only. Once 
spherical symmetry is broken, most radiation schemes be- 
come quite difficult to implement, given the large number 
of phase space degrees of freedom that need to be tracked 
for the radiation field. 

In this paper, we formulate the radiation transport 
equations in the framework of our 3 + 1 GRMHD scheme, 
which operates without any restrictions regarding the 
spatial symmetry of the system. However, our imple- 
mentation focuses on the optically thick limit for the ra- 
diation field, which simplifies the analysis by allowing us 
to assume that the radiation field in the comoving frame 
of the fiuid is nearly isotropic. Our emphasis is geared to 
treating systems in which the radiation has a strong dy- 
namical influence on the matter flow and, in some cases, 
on the spacetime geometry itself. We are less concerned 
in this initial treatment with the radiation that escapes 
from the matter surface, or with the radiation spectrum 
measured by a distant observer. It is in the interiors of 
collapsing stars, neutron stars in merging compact bina- 
ries, and dense accretion disks orbiting black holes where 
the dynamical influence of the radiation field is likely to 
play its most significant role. In these interior regions 
the optically thick assumption should be quite reliable 
in many cases. In the implementation presented here we 
also adopt a grey -body opacity law, which, though simple, 
suffices to illustrate our method. However, the formalism 
makes no assumptions regarding the spatial symmetry of 
the matter source, radiation field, or spacetime. 

We present two sets of tests to check our new radia- 
tion GRMHD code. The first set of tests involves radi- 
ation shocks and nonlinear waves propagating in a fixed 
Minkowski spacetime. The second set of tests is the 
"thermal Oppenheimer-Snyder collapse" problem origi- 



nally proposed and solved by Shapiro [29|,|30|, wherein ra- 
diation propagates in a spherical spacetime that, though 
simple, is highly dynamical and characterized by a strong 
gravitational field (i.e. one in which a black hole forms). 
In both sets of tests, we compare our numerical results 
with analytic solutions and perform convergence tests. 

The structure of the paper is as follows: In Sec. [Hi we 
formulate the system of coupled Einstein-Maxwell-MHD- 
radiation equations in 3+1 form, with the Maxwell, MHD 
and radiation equations written in conservative form. In 
Sec. IIII[ we describe techniques for evolving this system 
of equations. In Sec. |TVl we present the new code tests 
and their results. Finally, we summarize our results in 
Sec.lVl 



II. FORMALISM 

Throughout this paper, Latin indices denote spatial 
components (1-3) and Greek indices denote spacetime 
components (0-3). We adopt geometrized units, so that 
G = c= 1. 



A. Evolution of gravitational fields 

We write the spacetime metric in the standard 3+1 
form: 



ds" 



a'dt^ + 7y (dx' + (3'dt){dx^ +(3'dt), (1) 



where a, /3\ and 7^ are the lapse, shift, and spatial met- 
ric, respectively. The extrinsic curvature Kij is defined 
by 



(2) 



where £^3 is the Lie derivative with respect to The 
evolution of 7ij and Kij is governed by the Einstein equa- 
tion G^i, = SttT^i,, where G^i, is the Einstein tensor and 
Tfj^iy is the stress-energy tensor. 

We evolve 7^ and Kij using the BSSN formulation [l^, 
1 191. The fundamental variables for BSSN evolution are 



ln[det(7y)] 



K = 
Ai-j = 

r = 



1 

12 



(3) 

(4) 
(5) 

(6) 

(7) 



The Einstein equation Gfn, = SttT^^ gives rise to the 
evolution equations and constraint equations for these 
fields, which are summarized in [Tgj . In this paper, we 
use the same field evolution equations as Eqs. (11)-(15) 
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of [41|: 

{dt - Cp)^ 



'2a A, 
1 



1 



(8) 
(9) 

{dt-Cfj)K = -fWjD,a + ^aK'' (10) 

{dt - C0)A,, - e-''^{-D,D,a + a{R,, - 8^5,,))^^ 
+a{KA,j - 2AaA',), (11) 

and 



dtr = d,{2aA'^ + Cpf^) 



(12) 



2~ 



r(3\j+(3^T\j-2A'^d,a 



where D denotes covariant derivative operator associated 
with 7y , and TF denotes the trace-free part of a tensor. 
The constraint equations, expressed in terms of the BSSN 
variables, are 



= H = fWiD^e'^ - —R 



(13) 



=50 



-A,,A'^ - —K^ + 2TTe"^p, 



z^'I'D'K -8ne^'f'S\{U) 



where D denotes covariant derivative operator associated 
with . The matter-energy source terms are given by 



Si = -JiaTlpT"^ 

s = Y's.j . 



(15) 



Here n" = (a^^, — a^^/3*) is the time-hke unit vector 
normal to the t = constant time slices. In this paper 
T"^ contains three components: 



(hydro) (cm) 



(16) 



where T,?^ . s , T?'^ ^ and R""^ are the stress-energy ten- 

(nydroj ' (cm) ^^-^ 

sor for the hydrodynamic matter field, (large scale) elec- 
trodynamic field and the radiation field, respectively. 
Hence all components here contribute to the BSSN source 
terms in Eq. 

In order to evolve the 3-1-1 Einstein equations forward 
in time, one must choose lapse a and shift /3* functions, 
which specify how the spacetime is foliated. The lapse 
and shift must be chosen in such a way that the total 



system of evolution equations is stable. In the past few 
years, we have experimented with several gauge condi- 
tions. We find that, in general, the most useful gauge 
choices are the hyperbolic driver conditions fZ?, 43] , and 
the puncture gauge conditions (see e.g. [44j 45])- In this 
paper, we use the hyperbolic driver conditions as in j43j 
when evolving a dynamical spacetime: 



dtOi = aA 



dtA 



-ai {adtK + a2A -f a^^e-'^'I'aK) . (17) 



d'^P' = biadtV - b2dtl3' 



(18) 



where oi, 02, 03, 5i, and 62 are freely specifiable con- 
stants. 



B. Evolution of radiation fields 

1. Radiation fields 

The equations governing the dynamics of the radiation 
can be expressed as 



R 



;/3 



-G° 



(19) 



where i?"'^ is the radiation stress-energy tensor, and G" 
is the radiation four-force density which describes the 
interaction of the matter with the radiation J^] . The 
radiation stress-energy tensor R"/^ is defined as 



R 



a/3 



dvdVlhN'^N'^ 



(20) 



where v is the frequency, = /(x"; A^*, v) is the specific 
intensity of radiation at moving in the direction TV" = 
jf^ jhv^ is the photon 4- momentum, h is the Planck 
constant, and dfi is the differential solid angle. Here i/, 
li, and d£l are all measured in the local Lorentz frame 
of a fiducial observer with 4- velocity w^g^j), i.e. hv = 
~Pau"f^^y The integral is evaluated over all frequency 
and solid angles. 

We now choose our fiducial observer to be comoving 
with the fluid. In the comoving frame of the fluid the 
radiation stress-energy tensor R"^ takes the form 



R 



a/3 



£; pv pz 

jpx j^xx J^xy J^xz 

py 'pyi 'pyy 'pyz 

pz -pzx -pzy -pzz 



(21) 



where 



(22) 



E = j dvdVLl^ 
is the comoving radiation energy density, 

F'^ = / dvdni^N' (23) 



4 



is the comoving radiation flux, and 



dvdQ.I„N'N^ 



(24) 



is the comoving radiation stress tensor. 

We are interested in the optically thick regime, in 
which the radiation is very nearly isotropic in the co- 
moving frame of the fluid. In the limit of strict isotropy, 
independent of the propagation direction , the inten- 
sity is = I{x°'\ v). Using this fact and the expression 
of N°' in the comoving frame 

iV" = (1, TV*") = (1, sinOcofiLp, sin6lsinv3, cos6'), (25) 

one can show that = and V''^ = \5^'^E = S^^V, where 
■p is the radiation pressure, 6 is the polar angle measured 
from the ^-axis, and is the azimuthal angle (i.e. tantp = 
/N^) . Henceforth, we include the effect of a small 
anisotropy by allowing a small non-zero radiation flux 
F', but we retain the closure relation V = E/3. That is, 
we adopt an Eddington factor equal to 1/3. 

The radiation stress-energy tensor R"^ can be written 
in covariant form as 



R 



(26) 



where u" is the fluid 4-vclocity. This expression reduces 
to the same form as Eq. (PT|) in the comoving frame. Here 
we have introduced the projection tensor, h"^ , defined as 

and the radiation flux four-vector defined as 

F" = /i";3 J dvdVLl^N^. (28) 

Note that with this definition, the fiux satisfies 

F"ua = 0. (29) 
Following [s^l, the radiation four- force density is given 

by 



(30) 



where Xi^ — + xt is the total opacity (the superscript 
a and s denote the absorption and scattering opacities 
respectively) and rj^ — r]'^ + rjf, is the total emissivity. By 
assuming isotropic and coherent scattering, and that the 
thermal emissivity 77° and absorption coefficient x° are 
related by Kirchhoff's law = xtB,^, we can write, in 
the ffuid comoving frame. 



G" 



dvdnixlh J dvdnxlilu - B,) , 

dvd^lixl + xDhN', (31) 



where Bi, is the intensity in thermal equilibrium (e.g. 
the Planck function for photons, the analogous Fermi- 
Dirac function for neutrinos, etc.). We further assume 



a grey-body form for all opacities, Xv = '^Po, where k. 
is a frequency independent opacity, and pg is the rest 
mass-energy density. Then we may write [30| 



(32) 



It is straightforward to express G" in covariant form as 

= poK^iE - AttB)u°' + poiK" + k')F" . (33) 

Note that the frequency integrated equilibrium intensity 
B{T) can be written as 



47rB = aRT'^ 



(34) 



where T is the temperature of the fluid, and an is a con- 
stant depending on the type of radiation: for thermal 
photons it equals the usual radiation constant a; for each 
flavor of nondegenerate thermal neutrino or antineutrino 
(chemical potentials = 0) it is (7/16)a; lumping the con- 
tributions of all neutrinos and antineutrinos together, it 
is (7A/'iy/8)a, where AC is the number of neutrino flavors 
which contribute to thermal radiation. 

We emphasize here that our method allows for situa- 
tions in which the gas may be out of thermal equilib- 
rium with the radiation [E ^ AnB). Our formalism 
is equivalent to keeping the first two radiative moment 
equations, and using an Eddington factor to close the 
set. Our choice of V — 1/3E serves as the necessary 
closure relation for these equations. We demonstrate 
in Appendix lA II that our formalism, while more gen- 
eral, reduces to the relativistic diffusion approximation 
in a simplifying limit. In this limit, the diffusion ap- 
proximation transforms the radiation moment evolution 
equations from a hyperbolic to a parabolic (i.e. diffu- 
sion) form, which does not have the same causal struc- 
ture as the original system of equations. Moreover, the 
parabolic form is not suitable for implementing the con- 
servative HRSC scheme used to integrate the combined 
MHD-radiation equations (see Sec. IIIip . In any case, we 
do not adopt the diffusion approximation here, but treat 
the original set without simplification. 



2. Radiation evolution 

We can decompose the radiation evolution equations 
given by (jl9[) in a manner analogous to the way we de- 



corapose the MHD evolution equations (e.g., see Sec. IIC 
in TT'I and Sec. Ill Dl below). The resulting equations are 
therefore cast in conservative form, as are the MHD evo- 
lution equations. Taking the scalar product of Eq. (fT9)) 
with Ua on both sides gives the energy equation 



(35) 



where the radiation energy density variable r is deflned 
as 



f = ia^^)R' 



00 
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= Vl{auy^E + 2^a'u''F''~^^E , (36) 
and the source term s is 

+ . (37) 

Here 7 — e^^"^ denotes the determinant of the spatial 
metric 7ij . The spatial components of Eq. (|19p give the 
momentum equation, 

dtS^^d,{a^iRh) = aV7 Qi^^^'sa/s,. - ,(38) 

where the radiation momentum density variable is de- 
fined as 



is the magnetic field measured by an observer comoving 
with the fiuid. The stress-energy tensor associated with 
the perfect fiuid can be expressed as 



T; 



(hyd: 



fJ.V 



(44) 



where po is the (baryon) rest-mass density, P is matter 
pressure, h ^ 1 + e -f P/ po is the specific enthalpy, and e 
is the specific internal energy density of the matter. For 
brevity, we denote 



T, 



mhd) 



= t; 



Q/3 



(hydro) 



-t; 



Q/3 



(cm) 



(45) 



Thus, we see that the conservation of the total stress- 
energy tensor can be written as 



T 



■a(3 



t: 



Q/3 



(mhd) 



R 



a(3 



= 



(46) 



Si ^ ay^R°,i, 

= (^Eu°u, + F''u, + F,u''] . (39) 



C. Evolution of large-scale electromagnetic fields 

The evolution equation for the electromagnetic field 
in a perfectly conducting MHD fiuid {F'^'^Ui, = 0) can 
be obtained in conservative form by taking the dual of 
Maxwell's equation F[^^ x] = 0- One finds 



1 



{a^ *F^n = , 



(40) 



where F"^ is the Faraday tensor, and *F°'^ = 
e"^^'^Ff^i^/2 is its dual. Using the fact that the magnetic 
field as measured by a normal observer is given by 
= n^*F^^^, the time component of Eq. (|^ gives the 
no-monopole constraint djB^ = 0, where = ^ B^ . 
The spatial components of Eq. (|40|) give the magnetic 
induction equation, which can be written as 



dtB^ + d,j (v^B' -v^B^^^O 
where = /vP . 



(41) 



D. Evolution of the MHD field 



In the MHD limit, T^^^^-^ can be expressed as 



(cm) 2 



where h^^ = B^^^ / v47r and where 



^(«) —^v^ 



(42) 



(43) 



This can be combined with (flQl) to give 



T, 



ahd);/3 



(47) 



Additionally, we have the continuity equation expressing 
baryon number conservation. 



{pou'')-y = 



(48) 



Rewriting Eqs. (|47)) and ([48]) in conservative form gives 
(cf. Section HC in [l7|) 



dtp* + dj{p*v') = , 



dtS, -f 9j(av^T('^i^^)J = ^aN/7?^(m1id)5a/3,» + a^/lG^ , 

(50) 

dtT + a.(a2y^T°J,hd) - P*v') = s + c^VtG" , (51) 



(49) 



where the MHD evolution variables are 



(52) 



= ay/jpou° 

Si = \/7"-A''^(mhd)i 
= "V7?'(mhd)i 

= {p*h + au°y/jb'^)ui — ay/^b°bi , (53) 

T - V7'^M"^'?'(mhd) -P* 

= a'VlT^^i.d)-p* , (54) 

and the source term s is 
s = -aV7T('^rhd)^-% 



= [(7^(mhd)/3^/3^ + 2T('::,hd)/3^ + T^^nhd))^^. 



(^(mhd)^' + ^(mhd))t^i 



(55) 



Note that these evolution variables are very similar to 
those in [l^]. The only difference is that there are new 
radiative source terms Gi and G° in the momentum and 
energy equations (|5D|) and ([5T|l . respectively. 
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To complete the system of equations, it remains only 
to specify the equation of state (EOS) of the fluid. In 
this paper, we adopt a F-law EOS, 



P= (r-l)poe, 



(56) 



where F is the adiabatic gas constant. We choose a F-law 
EOS because it simplifies some of the calculations, it is 
applicable to many cases of interest, and it is a standard 
choice for demonstrating new computational techniques 
in the numerical relativity literature. Also, the analytic 
solutions we are going to use as code tests also use this 
EOS. Nevertheless, all evolution equations derived in this 
section apply for any equation of state, and generalization 
to a more realistic EOS is straightforward. In fact, our 
code is currently capable of handling the general class of 
EOSs of the form P = P{po,e). 

The fluid temperature T is required in the radiation 
force density term (Eq. ([34]) ). In this paper, we com- 
pute it by using the ideal gas law P = nksT = pokBT/m, 
where n is the baryon number density, m = po/n is the 
mean mass of the baryons in the fluid, and fcs is Boltz- 
mann's constant. Hereafter we set ks = 

We point out that the fluid flow is nonadiabatic in gen- 
eral. In particular, there is energy exchange between the 
matter and radiation fields. Also, shocks may be present 
in some applications. 



E. Summary of equations 

To reiterate, the system of coupled Einstein-radiation- 
Maxwell-MHD equations we consider are the BSSN equa- 
tions ((8))~ p^ . the radiation transport equations p5|) and 
([38|) . the magnetic induction equation (|41]) . and the MHD 
equations (|49 |) --([5T |) . In Appendix El we demonstrate 
that our equations reduce to the more familiar Newto- 
nian form in the weak-field, slow-velocity limit. The evo- 
lution variables are 0, 7^, Aij, F% f. Si, B^, p*, Si 
and f. These variables are not completely independent: 



the BSSN variables 0, 7^^, K, Aij, and F' have to sat- 



isfy the Hamiltonian constraint (jl3p and the momentum 
constraint (fT4|) : the magnetic field variables have to 
satisfy the no-monopole constraint 9^5* = 0. 
The total stress-energy tensor T'"' is given by 



(hydro) (cm) 



poh + + -Ej uf^u" 



The BSSN matter-energy source terms [Eq. 
expressed as 



{au°f ( p^h + b^ + -E 



(57) 



ISllI can be 



-E 



S, = 



+2a^u'>F° - {ab"f (58) 
auo (^poh + fe^ + ^E^ u, + aF"ui + au°F^ 



-ab\ (59) 
pnh + b^ + ^Ej u,uj +ip+ + ^Ej -f,j 



-FiUj + FjUi - bibj 



III. IMPLEMENTATION 



(60) 



We use a cell-centered Cartesian grid in our three- 
dimensional simulations. Sometimes, symmetries can be 
invoked to reduce the integration domain. For octant 
symmetric systems, we evolve only the upper octant; for 
equatorially symmetric systems, we evolve only the up- 
per half-plane. For axisymmetric systems, we evolve only 
the x-z plane (a 2-f 1 dimensional problem). In axisym- 
metric evolutions, we adopt the Cartoon method [13] for 
evolving the BSSN equations, and use a cylindrical grid 
for evolving the induction, MHD, and radiation equa- 
tions 

Our code uses the Cactus parallelization framework 
[4^, with the time-stepping algorithm based on the 
MoL, or method of lines, thorn. In the metric evolu- 
tion (BSSN sector), spatial derivatives can be calcu- 
lated using second-order or fourth-order finite differenc- 
ing schemes. The Cactus MoL thorn allows us to switch 
to a higher-order time-stepping scheme easily. Higher or- 
der schemes are very useful for evolving spacetimes con- 
taining black holes using the moving puncture techniques 
(see, e.g., [13, Iks.]). However, we do not treat punc- 
ture black holes here and we are currently using a HRSC 
scheme which is at most second-order accurate to evolve 
the Maxwell, radiation, and MHD equations. Hence we 
use second-order finite differencing scheme in the BSSN 
sector and (second-order) iterated Crank-Nicholson time- 
stepping in our calculations. 

Our technique for metric evolution is described in our 
earlier papers , so we focus here on our MHD, 

induction and radiation algorithms. The goal of this part 
of the numerical evolution is to determine the fundamen- 
tal "primitive" variables P = {po, P,v^ , B^ , E, F^) at fu- 
ture times, given initial values of P. The evolution equa- 
tions dini), (IMl), 611), (HSl), (ISO]) and ((51]) are written in 
conservative form: 



9tU + V • F = S 



(61) 



where the evolution variables U(P) = (p* ,f ,S'i,i3',f,5i), 
the fiuxes F(P) and the sources S(P) are not explicit 
functions of derivatives of the primitive variables, al- 
though they are explicit functions of the metric and its 
derivatives. As mentioned above, we evolve Eq. (j6ip us- 
ing the iterated Crank-Nicholson scheme. This scheme 
is second order in time and will be stable if At < 
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niin(Aa;*)/cmax, where in our case Cmax is the speed of 
light. For each Crank-Nicholson substep, we first update 
the gravitational field variables (the BSSN variables). We 
then update the electromagnetic fields by integrating 
the induction equation. Next, the MHD variables (p*, f , 
and Si) are updated. Then we update the radiation vari- 
ables (f, and Si). Finally, we use these updated values 
to recover the primitive variables on the new timestep. 
Below, we briefly summarize some of the important tech- 
niques we utilize during the evolution. 



as described in Appendix |B] To obtain us/ki in the 
grid frame, we use the dispersion relation (|B4p and 
substitute the ojcm = —k^j,u^^, and fc^jj, = K^K^, where 
= {gp,i, + u^u^)k'^ . Wave speeds in the y-direction 
and z-direction are found analogously. 



C. Recovery of primitive variables 



A. Reconstruction step 

We implement an approximate Riemann solver to han- 
dle the advection in Eq. (pT|) . For simplicity, we con- 
sider the one-dimensional case here. The generalization 
to multi-dimension is straightforward. The first step in 
calculating this flux is to compute Pl = Pi+i/2-£ and 
P/j = Pi+i/2-i-ej the primitive variables to the left 
and right of the grid cell interface. As in we use the 
Monotonized central (MC) scheme [5l| to compute the 
primitive variables at the cell interface. This scheme is 
second-order accurate at most points when the data are 
smooth, but becomes first-order accurate across a dis- 
continuity (e.g. shock). (See [l7| for other reconstruction 
methods in our MHD code.) 



B. Riemann solver step 

Next, we take the reconstructed data as initial data for 
a piecewise constant Riemann problem, with P — P^ on 
the left of the interface, and P = P^ on the right of the 
interface. The net flux at the cell interface is given by 
the solution to this Riemann problem. 

We use the HLL (Harten, Lax, and van Leer) approxi- 
mate Riemann solver [52] . Our implementation has been 
described in [l7[. To summarize, HLL fluxes are given 
by 



/'i+l/2 

Here 



--min'.'max 



(ur - Ul) 



+ c„ 



Cmax = max(0, C+l) 
Cmin = -min(0, c_fl,,c_i) 



(62) 



(63) 



where c+ is the maximum right-going wave speed and 
c_ is the maximum left-going wave speed. We obtain c± 
by solving the dispersion relation for waves with wave 
vectors of the form 



fc^ = (-a;,fci,0,0) 



(64) 



The wave speed is simply the phase speed w/fci. We 
find the dispersion relation in the comoving frame of the 
fluid (denoted by subscript cm), and hence Wcm/fccm, 



Having computed U at the new timestep, we must use 
these values to recover P, the primitive variables on the 
new time level. We can recover the hydrodynamics prim- 
itive variables po. Pi v'^ from the MHD evolution variables 
p*, f , Si numerically, as described in Section HI C in p^j . 
Once the fluid velocity is found, the radiation primitive 
variables E and F"^ can be computed from the radiation 
evolution variables f and Si analytically using Eqs. (|36p 
and ([55]) . We solve the following set of two coupled linear 
equations to recover E and F^: 



= V7 



-(au°)2- i ) £; + 2a'u"F^ 



(65) 



- au°f + {u^P' -t- u')S^ = -a^{Eu° + F°) (66) 

The first one is just Eq. ((36)) . while the second one is 
obtained by using u^F^ = to eliminate Fi in Eq. (|39|1. 
After solving for E and we compute by 

F' ^ - -Eu° (v' + f3') - 2F°I3' ~ F%' , (67) 

a^u^' 3 ^ ' 

which is derived by raising the index of Fi in Eq. (j39p . 

D. Constrained Transport 

The Maxwell equation demands that the magnetic 
fields satisfy the no-monopole constraint diB'' ~ 0. 
Unphysical behavior may arise if this constraint is vio- 
lated. Thus, "constrained transport schemes" have been 
designed to evolve the induction equation while main- 
taining diB^ = to roundoff precision ^5^. We use the 
flux-interpolated constrained transport (flux-CT) scheme 
introduced by Toth [20] and used by Gammie et al [s^ ]. 
This scheme involves replacing the induction equation 
flux computed at each point with linear combinations 
of the fluxes computed at that point and neighboring 
points. The combination assures both that second-order 
accuracy is maintained, and a particular finite-difference 
representation of diB^ — is enforced to machine prcci- 
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E. Low-density regions and boundary conditions 

1. Low- density regions 

If vacuum exists anywhere in our computational do- 
main, the MHD approximation will not apply in this 
region, and we will have to solve the vacuum Maxwell 
equations there (see e.g. [111). In addition, the optically 
thick assumption on the radiation field in Sec. Ill Bl also 
breaks down in sufficiently low density regions. In many 
astrophysical scenarios, however, a sufficiently dense, ion- 
ized plasma will exist outside the stars or disks, whereby 
MHD will remain valid in its force- free limit. A sim- 
ilar situation may arise for the radiation field, where 
the ambient gas in our computational domain may be 
sufficiently dense to maintain an optical depth above 
unity. However, in some applications we may need to 
take into account the transition from optically thick to 
optically thin limits in the low-density regions, depend- 
ing on the magnitude of the opacity. A precise treatment 
of this problem requires solving the full Boltzmann ra- 
diative transfer equation (see e.g. [13, [H, [s^, [Ell; See 
also [H, [ill for approximation schemes.) In this pa- 
per, however, we avoid this issue. For the code tests 
that do not have low density regions (Sec. lIVX)) . no spe- 
cial treatment is required. We do, however, present a 
test involving Oppenheimer-Snyder collapse (Sec. lIVBp 
where there is a vacuum outside the star. As in many 
hydrodynamics simulations in astrophysics, we impose a 
low-density "atmosphere" outside the star to facilitate 
the integration of hydrodynamics equations. It turns out 
that our atmosphere scheme suffices to mimic the cor- 
rect ( "zero temperature" ) radiation boundary conditions 
that we wish to impose at the surface of the star (see 
Sec. UVBl) . 

In the low-density regions near the surface of the 
star, we sometimes encounter problems when recover- 
ing the primitive variables; in particular, the equations 
U = U(P) occasionally have no physical solution. Usu- 
ally, unphysical U are those values corresponding to neg- 
ative pressure. As in [l3|, we apply a fix at these points, 
first suggested by Font et al [6i|. In the system of equa- 
tions to be solved, we replace Eq. with 
the adiabatic relation P = kpq, where k is set equal to 
its initial value. This substitution guarantees a positive 
pressure. Typically, these low density regions have lit- 
tle influence on the dynamical evolution of the system, 
which is the principal target of our current investigations. 

2. Boundary conditions 

For the code test in Sec. IIV A| we evolve one- 
dimensional, radiation-hydrodynamics equations in a 
fixed, Minkowski spacetime. We impose the "copy" 
boundary condition on all the evolution variables, i.e. 
variables at the boundaries are copied from the closest 
grid point. 



For the Oppenheimer-Snyder code test in Sec. IIVBI 
we evolve the system of coupled Einstein-radiation- 
hydrodynamics equations. In this case, we employ Som- 
merfeld outgoing wave boundary conditions for all BSSN 
and gauge variables: 

/(r,t) = ^l^/(r-Ar,i-AT), (68) 
r 

where AT is the timestep and Ar = ae~^'*AT. For the 
radiation hydrodynamics, we impose the outflow bound- 
ary condition on the primitive variables p, P, w*, E, and 
(i.e., the variables are copied along the grid directions 
with the condition that the velocities be positive or zero 
in the outer grid zones). We note that the radiation in 
this test is initially confined inside the star, but escapes 
from the stellar surface during the evolution. In the end 
of the simulation, the total emitted radiation remains 
small and the dynamics of the system is insensitive to 
the boundary condition employed. 

IV. CODE TESTS 

Our GRMHD code has previously been thoroughly 
tested by maintaining stable rotating stars in stationary 
equilibrium, by reproducing Oppenheimer-Snyder col- 
lapse to a black hole, and by reproducing analytic solu- 
tions involving MHD shocks, nonlinear MHD wave prop- 
agation, magnetized Bondi accretion, and MHD waves 
induced by linear gravitational waves [l3|. It has also 
been compared with the GRMHD code of Shibata and 
Sekiguchi |62] by performing simulations of the evolu- 
tion of magnetized hypermassive neutron stars |3j8| , and 
of magnetorotational collapse of stellar cores [2l|. We 
obtain good agreement between these two independent 
codes. Our code has also been used to study the evolu- 
tion of BHBH and BHNS binaries [23| , and the evolution 
of relativistic hydrodynamic matter in the presence of 
puncture black holes [63] . Here we restrict our attention 
to testing the new radiation-hydrodynamics sector, set- 
ting large-scale magnetic fields to zero. We also choose 
the grey body absorption opacity k° to be a constant and 
set the scattering opacity to zero. 

A. Minkowski radiation-hydrodynamics tests 

We present here a series of tests of nonlinear radiation- 
hydrodynamic waves in Minkowski spacetime with pla- 
nar symmetry. These tests are summarized in Table ID 
For our initial data, we generate semi- analytic, stationary 
configurations using the method outlined in Appendix [Cl 
To test the ability of our code to handle shocks and waves 
moving across the grid, we boost the stationary solutions 
derived in Appendix [C] for tests 2-4. In each case, our 
computational domain is a; G (—20,20). We choose the 
opacities in each case to ensure that the grid boundaries 
at a; = ±20 reside in the asymptotic region where all 
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TABLE L Initial states for one-dimensional tests. 

Test f 7? Left state'' Right State'' tenai tsc" 

T 5/3 Oi pa = 1.0 po = 2.4 5000 2000 

(/i = 0.0)=' P = 3.0 X 10"^ P = 1.61 X 10-'' 

= 0.015 u'' = 6.25 X 10"^ 

£ = 1.0 X 10"^ E = 2.51 X 10"^ 

"2 5/3 02 po = 1.0 po = 3.11 100 80 

= 0.1)'' P = 4.0 X 10-=* P = 0.04512 

= 0.25 u"" = 0.0804 

g = 2.0 X 10"° E = 3.46 X 10"^ 

"3 2 03 po = 1-0 po = 8.0 20 20 

(p = 0.8)'' P = 60.0 P = 2.34 X 10^ 

u"" = 10.0 = 1.25 

E = 2.0 E = 1.14 X 10^ 

4 5/3 0.08 po = 1.0 Po = 3.65 100 90 

(p^O.l)'' P = 6.0 X 10"^ P = 3.59 X 10"^ 

u"" = 0.69 u'' = 0.189 

£ = 0.18 i; = 1.30 



'^ /i is the speed at which the wave travels. Traveling wave solutions are obtained by 
boosting the stationary solutions in Appendix [Cl to speed p. 
^ tac is the approximate time it takes for a wave traveling at the sound speed to 
propagate from the center of the grid to the right boundary. 

Values refer to asymptotic regions. We solve ODEs to determine the exact 
solution in the transition region (see Appendix [C|) . 



hydrodynamic and radiation quantities approach their 
asymptotic values, and that the total optical depth across 
the grid is r ^ 10 (see Appendix [Cjl . 

We evolve the system with a timestep At = Ax for test 
1 and test 2, and At = O.lAx for test 3 and test 4. We 
use resolutions ranging from Ax = 0.0125 to Ax — 0.1 
in order to perform convergence tests. To demonstrate 
convergence, we consider a grid function g with error 
Sg = g — g'^^'^'^^. We calculate the LI norm of Sg (the 
"average" of Sg) by summing over every grid point i: 

N 

Ll{6g)^AxJ2\g,-g^^'''^\x.)\ , (69) 
1=1 

where N oc I /Ax is the number of grid points. We find 
that for our continuous configurations (tests 3 and 4), 
we achieve second order convergence. Because our shock 
capturing scheme becomes 1st order when discontinuities 
are present, we achieve the expected first order conver- 
gence for our discontinuous configurations (tests 1 and 
2). 

The initial configurations listed in Table U are chosen 
to test our code in a variety of regimes, including gas- 
pressure dominated, radiation-pressure dominated, New- 
tonian, relativistic, continuous, and discontinuous matter 
and radiation profiles. In each test, the equation of state 
of the gas is given by a F-law EOS. We choose F = 5/3 
for each test except our highly-relativistic case, in which 
we choose F = 2. This latter choice is adopted because 
the sound speed (cg = ^TP/poh) for a F-law EOS is 
limited by Cg < \/r — 1. Highly-relativistic sound speeds 



{cs 1) can only be achieved for F > 2. Below, we 
provide a brief description of each test. 

1. Nonrelativistic strong shock. For this test, we set 
up a strong, gas-pressure dominated, Newtonian 
(■"max = 0.015 ^ 1) shock propagating into a cold 
gas. We have chosen to simulate this scenario be- 
cause it can be compared to the analytic solution 
for a subcritical radiating shock first derive d by 
Zel'dovich and Raizer [64] and summarized in [46|. 
We find very good agreement with this analytic re- 
sult (see Fig. [1]). We note that the radiative shock 
junction conditions (see Appendix [C]) require that 
R^^ and be continuous at the shock front, even 
though E and F'^ are, in general, discontinuous at 
the shock. In the Newtonian limit, however, the 
continuity of and R'^^ is equivalent to the con- 
tinuity of E and F^. 

2. Mildly -relativistic strong shock. In this test, we set 
up a mildly-relativistic (uJiax — 0.25), gas-pressure 
dominated shock. In this case, we see that E and 
F^ no longer appear continuous. We boost this 
shock so that the shock speed is fi — 0.1. We find 
that the discontinuity is able to retain its shape 
very well as the shock travels and matches very 
well with the analytic solution (see Fig. [J). 

3. Highly-relativistic wave. In this test, we simulate a 
highly-relativistic (w^ax = 10), gas-pressure domi- 
nated configuration in which all quantities are con- 
tinuous, but asymptote to different values on either 
side of the computational domain. We boost this 
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FIG. 1: Profiles of po, P, v'^ , E, and F"" &t t = 5000 for test 
1. In this test, the shock front remains stationary. Solid dots 
denote data from numerical simulations with resolution Aa; = 
0.0125. Solid lines denote the exact solutions (Appendix [Cj. 



FIG. 2: Profiles of po, P, v"" , E, and F"" at t = 100 for 
test 2. In this test, the shock front moves with velocity /i = 
0.1. Solid dots denote data from numerical simulations with 
resolution Ax = 0.0125. Solid lines denote the exact solutions 
(Appendix [C}. 



configuration so that it travels across the grid with 
velocity /i = 0.8. The numerical results agree very 
well with semi- analytic solution (see Fig. [3]) . Fig- 
ure |4] shows the LI norms of the errors in F'^, E, 
v'^, P and pq a.t t ~ ifinai ~ 20. We find that all 
errors converge to zero at second order in Ax. 

4. Radiation-pressure dominated, mildly-relativistic 
wave. In this test, we study the performance of our 
code in the radiation-pressure dominated {P <^V), 
mildly-relativistic [u^^^ = 0.69) regime. We boost 
this configuration so that it travels across the grid 
with velocity = 0.1. The numerical results again 
agree with the semi-analytic solution (see Fig. [5]) . 



B. Dynamical Spacetime test: Thermal 
Oppenheimer-Snyder Collapse 

The collapse from rest of a homogeneous dust ball 
[P ~ 0) in general relativity can be described by the an- 
alytic Oppenheimer-Snyder solution [65*1 . The collapse 
results in the formation of a Schwarzschild black hole. 
The evolution of thermal radiation within the dust ball 
has been considered by Shapiro in [2^, [s^ . In both pa- 
pers, the radiation is assumed to be a small perturbation, 
so that the dynamics are unaffected by the presence of 
radiation, and the matter and metric profiles can still be 
described by the Oppenheimer-Snyder solution. The first 
paper employs the relativistic thermal diffusion approx- 
imation for the radiation, and derives analytic solutions 



for both the Newtonian and general relativistic cases. 
In the second paper this approximation is removed and 
replaced by solving the exact radiative transfer (Boltz- 
mann) equation for the intensity, coupled to the radia- 
tion moment equations for the radiation flux and energy 
density. It is found that the results obtained by solv- 
ing the Boltzmann transport equation agree very well 
with the analytic solutions assuming diffusion approxi- 
mations, provided the optical depth of the star is suffi- 
ciently large (^ 1). Here we perform a numerical simu- 
lation of "thermal Oppenheimer-Snyder collapse" using 
our radiation GRMHD code, and compare our results to 
the analytic solution in the diffusion approximation limit 
given in (29l . [30| . For convenience, we summarize the 
analytic solution in Appendix [P] 

For our initial data, the areal radius of the star is set 
to be Ri = 3M, where M is the ADM mass of the star. 
We choose the initial profiles for all hydrodynamic and 
radiation quantities to be homogeneous throughout the 
star, in compliance with the analytic solution in Ap- 
pendix |D] The analytic solution assumes that (1) the 
matter and radiation pressure is small enough to be dy- 
namically unimportant (i.e. P/pa <C M/R and V / po — 
E/3po <C M/R), (2) radiation pressure dominates over 
gas pressure {V ^ P), (3) gas and radiation are in local 
thermal equilibrium (LTE), i.e. E = AttB — auT'^, and 
(4) the star is optically thick. To satisfy these conditions, 
we choose the following initial data: po = M/ {^irRf), 
P ^ IQ-^po, E = 10-^0, = 0, F' = 0. LTE 
is achieved in the initial data by fixing the constant 
flflm" = m^E/T^ = E{po/Py = IO^^M/HttRI). We 
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FIG. 3: Profiles of po, P, v"" , E, and f ^ at t = 20 for test 3. 
In this test, the shock front moves with velocity /i — 0.8. 
Solid dots denote data from numerical simulations with res- 
olution Ax = 0.0125. Solid lines denote the exact solutions 
(Appendix [C}. 



note that in our formalism the system is allowed to de- 
viate from the LTE during the evolution. However, we 
find that the system remains close to the LTE during the 
entire evolution and the numerical data agree well with 
the analytic solution (see below). We choose k"' so that 
r° = K°'pR = 50 initially. This guarantees that the star 
is optically thick initially. As the collapse proceeds, the 
optical depth increases as r cx so the star remains 

optically thick. 

We construct the initial data for the spatial metric 
by transforming the analytic Oppenheimer-Snyder met- 
ric from Friedmann to isotropic coordinates, following 
the procedure described in §31- We use the analytic so- 
lution only at t=0. The metric at later times is evolved, 
together with hydrodynamics and radiation. The lapse 
and shift are determined by the hyperbolic driver condi- 
tions (Eqs. P7)) and HI])). These are gauge conditions 
that have been widely used in stellar collapse calculations 
using the BSSN scheme. We choose oi = 0.75, &i — 0.15, 
a2 — b2 — 2M~^, 03 = 1. A smaller bi prevents "blowing 
out" of the coordinate system, a well-known effect [4l|, |67| 
which can spoil grid resolution in the center of the col- 
lapsing object. We perform our numerical simulation 
in axisymmetry, with 200^, 400^, 800^ and 1600^ grid 
points. We choose At = 0.1 Ax in these simulations. 
The outer boundary is placed at 4M in isotropic coordi- 
nates {Rout = 5.06M in areal radius). Note that we do 
not impose any special boundary condition at the stellar 
surface, in contrast to the zero temperature boundary 
condition {E = 0) used in the derivation of the ana- 
lytic solution [29j. The low density region outside the 



FIG. 4: LI norms of the errors in po, P, v"^ , E, and F"" for 
test 3 at i = 20. This log-log plot shows that the LI norms 
of the errors in all quantities are proportional to (Ao;)^, and 
are thus second-order convergent. 



star mimics this surface boundary condition, as the at- 
mosphere is made to be much colder than the interior of 
the star, and hence the thermal emission and build-up of 
radiation energy density in the atmosphere is negligible. 

The analytic solution given in Appendix|D]is expressed 
in Friedmann coordinates (i.e. Gaussian normal coordi- 
nates comoving with the fluid), which is equivalent to 
using the gauge conditions a = 1 and = (geodesic 
slicing and zero shift), which are different from the gauge 
conditions we adopt in our numerical simulations. In or- 
der to compare our numerical result to the analytic solu- 
tion, we perform a mapping between these two different 
gauges. This is achieved first by following a set of La- 
grangian fluid elements inside the star, and calculating 
the proper time and position of these elements by inte- 
grating the equations 



dr 
'dt 



1 



~dt 



(70) 



Next, we use the metric and the positions of the fluid el- 
ements to compute their areal radii r^. Finally, knowing 
the proper times t and areal radii Tg of the fluid elements, 
we use Eqs. (jD2|) and (ID3|) and rg = a(T) sin x to com- 
pute their Friedmann coordinates (r, x) ■ The mapping 
between these two gauges is thus established. The pair 
(tj, Xj) for each element j uniquely specifies the fluid and 
radiation parameters. 

Figures [H] and [7] show the profiles of po, P, E and F at 
different times during the collapse. Note that while the 
density remains spatially constant during the collapse in 
comoving Friedmann coordinates, this is not true in our 
gauge. We see that the numerical results agree very well 
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FIG. 5: Profiles of po, P, v"" , E, and F"" at t ^ 100 for 
test 4. In this test, the shock front moves with velocity ^ = 
0.1. Solid dots denote data from numerical simulations with 
resolution Ax — 0.0125. Solid lines denote the exact solutions 
(Appendix [C}. 
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FIG. 6: Profiles of the hydrodynamic quantities po and P in 
Schwarzschild (areal) radius at times t/M=Q, 6, 8 and 9.5 for 
thermal Oppenheimer- Snyder collapse. Solid lines represent 
numerical data with 1600^ grid points, and dashed lines show 
analytic solutions. The solid square on the a:-axis denotes the 
radius of the apparent horizon, which forms after the stellar 
surface passes through an areal radius of 2M. 



FIG. 7: Same as Fig.[6]but for the radiation quantities E and 
F inside the star. Note that F = G everywhere at t = 0. 

with the analytical solution, even after the apparent hori- 
zon appears at t = 6.78M and all of the stellar material 
is inside the horizon. 

We next perform a convergence test for the radiation 
quantities. We find that our numerical data converge to 
a solution slightly different from the analytic solution. 
This can be explained by the fact that the analytic solu- 
tion is strictly valid only in the perturbative limits that 
P/po ^ and P/V 0. While we set up our initial data 
to approximate these limits, the slight deviation from the 
analytic solution is still detectable with our resolutions. 
In the absence of radiation {E = 0, F'^ = 0), we have 
checked that the deviation in po between numerical and 
analytic values is reduced by a factor of 10 if we reduce 
the ratio P/ po by a factor of 10. In the presence of radia- 
tion, however, decreasing the ratios P/V and V/po arbi- 
trarily small makes the numerical simulations quite chal- 
lenging, as accurate evolution for the radiation quantities 
E and requires accurate evolution of the temperature 
T (X P/ Po, which in turn requires accurate determination 
of the pressure P. However, accurate computation of P 
in the limit P/po ^ is difficult. Since the evolution 
variables are dominated by the rest mass density po, in 
order to recover the tiny P accurately from them, the nu- 
merical evolution has to be very accurate. This requires 
very high resolution. Thus, we perform a convergence 
test in which we compare numerical solutions with small 
but finite P/V and P/ po for different resolutions, rather 
than comparing with the analytic solution. 

Figure [S] shows the result of the convergence test for E 
and F, with the differences scaled for second order con- 
vergence. We follow a Lagrangian fluid element halfway 
between the center and the surface of the star (in terms 
of areal radius), determine the radiation parameters at 
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FIG. 8: Convergence test for the radiation quantities E and 
F, computed at the Lagrangian point halfway between the 
center and the surface of the star (in terms of areal radius). 
Resolutions with 200^, 400^, 800^ and 1600^ grids are used 
here. The differences between lower and the highest solu- 
tion (1600^) are rescaled to demonstrate second-order con- 
vergence. 



the position of the element versus time, and subtract the 
numerical results from different resolutions. Since we use 
a HRSC scheme that is second-order accurate except at 
positions where discontinuities appear, e.g. at the sur- 
face of the star where the density falls abruptly, we ex- 
pect that the order of convergence depends on how much 
the physical quantity is affected by the discontinuity at 
the stellar surface. In principle, the first-order error will 
propagate everywhere inside the star, but its effect may 
be small (depending on the physical quantity under con- 
sideration) and may only be detectable at very high res- 
olution. Fig. [5] shows that E converges at second order, 
whereas F converges at less than second order but better 
than first order. This shows that F is more susceptible to 
propagation of the first order behavior at the surface of 
the star, which can be anticipated by looking at the shape 
of the profiles in Fig. [71 E drops abruptly before reach- 
ing the surface, while F increases monotonically up to the 
surface. We see that convergence of F deviates further 
from second order as the resolution is increased. This is 
consistent with the presence of a first order term with 
small coefficient due to the discontinuity at the stellar 
surface: F{A,t) = F^^^^t{t) + ci{t)A + C2{t)A'^ + 0{A^). 
Here F{A,t) is the value of F at time t evolved with 
a grid size A, Foxact(i) is the exact solution, and ci{t) 
andc2(i) are resolution-independent functions. The first- 
order term ci(t)A results from the discontinuity. We 
expect that Ci{t) <C C2{t) since we look at a point far 
away from the discontinuity. With a lower resolution. 



and hence a larger grid size A, the first-order term is not 
as significant relative to the second order term because of 
the small coefficient. Once we decrease A by increasing 
the resolution, the second-order term diminishes as A^, 
while the first order term shrinks as A only, making the 
first order term more conspicuous. Similar behavior for 
E is not evident since E drops to a very low value at the 
surface. 



V. CONCLUSIONS 

We have developed a code which can evolve the cou- 
pled Einstein-Maxwell-MHD-Radiation equations in 3 -|- 
1 dimensions. In the implementation presented here this 
code is able to model the behavior of magnetized, per- 
fectly conducting, radiating fiuids in dynamical space- 
times in which the characteristic length scales of the 
system are longer than the mean free path of the ra- 
diation and the opacity has a grey-body form. Our 
formalism allows us to evolve the radiation fields us- 
ing a HRSC scheme which is analogous to the method 
we use for the hydrodynamic fields. In this paper, we 
have tested the shock-capturing capabilities of our code 
by simulating both continuous and discontinuous one- 
dimensional radiating hydrodynamic waves. We have 
been successful in evolving highly relativistic, radiation- 
pressure dominated, gas-pressure dominated, and Newto- 
nian waves. We have treated both stationary waves and 
boosted waves that propagate across our computational 
grid in Minkowski spacetime. We have also confirmed 
our ability to accurately capture the behavior of radia- 
tion in a strong-field dynamical spacetime by simulating 
a thermal Oppenheimer-Snyder collapse. Our numerical 
results agree well with the analytic solutions. We per- 
form convergence tests on our test problems and find the 
expected order of convergence in all cases. 

We plan to use our radiation GRMHD code to study 
many interesting systems, and revisit some of the prob- 
lems we have considered before, such as core-collapse su- 
pernovae, accretion onto a black hole, merging NSNSs 
and BHNSs, etc. By taking into account the effect of ra- 
diation, we hope to gain more insights and provide some 
answers to questions that are relevant to observations. 
For example, we can calculate the radiation luminosity, 
and study the radiation feedback to the dynamics of the 
systems. 
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APPENDIX A: LIMITS 

1. Diffusion Approximation 

Using the form of R"^ given by Eq. ([26]) in the radia- 
tion moment equation R^^-p = —G" gives 



(Al) 



where V = E/3. In our formalism, we have already as- 
sumed near isotropy, which implies that F" /E <^ 1. We 
can make the further approximation that the radiation 
flux can be neglected in the left hand side of Eq. (|Aip 
to get 



;/3 



(A2) 



Operating on both sides of the above equation with the 
projection tensor h''a, using Eq. p3|) for G" and the fact 
that F^Ua — 0, we get 



3 ' 



Ea" + -E''°' 
4 



is-" 
3 



(A3) 



,0 



where k = k° -|- is the total opacity, a" = u'^-pu^ is 
the 4-acceleration, and where we have used the fact that 
h'^au" = and V = EjZ. Thus, we arrive at the expres- 
sion for the radiation flux in the diffusion approximation. 



1 



3 -I- K«)po 



h • ( \\1E + 3.E 



(A4) 



This is the relativistic diffusion equation relating the ra- 
diation flux i^" to the local radiation energy density E. 
If we further assume LTE, then E = qhT'^, in which case, 



^Xthh. ■ (VT + aT) , 



where 



A 



th — T, 



3 [k."- + k')pq 



(A5) 



(A6) 



This familiar result is in agreement with Eq. (3.2) of [29| . 
and with Eq. (2.5.28) of o8]. For a further discussion of 
the simplification in Eq. (jAip leading to the diffusion 
approximation, see [30j . 



2. Newtonian Limit 

It is instructive to consider the weak-field, slow- 
velocity (Newtonian) limit of the GR radiation hydro- 
dynamic equations, and show that our equations reduce 



to the familiar expressions of Newtonian radiation hydro- 
dynamics. For simplicity, we set all large-scale electro- 
magnetic fields to zero (i.e. T^^^^ = = B*). 



a. Continuity equation 

From Eq. (|48|). we have 

{poun-. = (AT) 

In Newtonian limit, the covariant derivative reduces to 
partial derivative (in Cartesian coordinates), and the 4- 
vector reduces to u" « (!,«*)• Hence the continuity 
equation reduces to the familiar expression: 

dtpo + d,ipov')^0 . (A8) 

b. Euler equation 
From Eq. we have 

at(aV7r°.) + dj{a^Th) = \a^T^^9c.p. ■ (A9) 

In Newtonian limit, the metric can be approximated by 

ds^ = -{l + 2<P)dt'^+{l-2<^>){dx^+dy^ + dz^) , (AlO) 

where $ <C 1 is the Newtonian gravitational potential. 
Keeping only the lowest order terms, we obtain 



(All) 



Using Eq. ([57|) and assuming po '3> P, po ^ E, and 
pov^ ^ F^, we obtain 

dtipov,) + d, [pov^v, + {P + V)5^A = -pom . (A12) 



Combining Eq. (|A12p with the continuity equation (jA8[) 
yields 

dtv, + v^djv, = -—d,{P + V)- 9.$ , (A13) 
po 

which is the familiar Newtonian Euler equation, allowing 
for gas plus radiation pressure, {P + 7^). 



c. Energy equation 

The energy equation in the Newtonian limit is derived 
by contracting with the equation ^ JT^'^ — 0. Using 
Eq. (|46)) and the continuity equation V^(/9ow^) = 0, we 
find, after some algebra, that 

w^V^(poe + F) + (poe + ^ + P + P)V„u'^ + V^P'' 
-t-F^a^ = . (A14) 
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where a^^ = u'^W^u^ is the 4-acceleration. We note that 
M^flp — and F'^Ufj^ — 0, which imphes that 



In Newtonian Hmit, 



F^a^ = F'a, (l + 0{v^)) 



Furthermore, ^ Fw' and ~ v^djVi, so that 



(A15) 



(A16) 



so we may neglect F'^a^ in the Newtonian hmit. Also, 
dtF° « dt{viF^) ~ dt{Ev^) <C 9tF so we may neglect 
this term as well. 

Hence in the Newtonian limit, we obtain 

dtipoe + E)+ v'd,{poe + E) + {poe + E + P + V)dy 
+diF' = . (A17) 



S^ 
f 

Rh 



E , 
F\ 



(A25) 
(A26) 
(A27) 
(A28) 



Inserting these into (jA23[) and (|A24[) . and dropping 
higher order terms, we obtain 

dtF, + d^r = -G„ (A29) 
dtE + djF^ = -G°, (A30) 

which agree with Eqs. (94.2) and (94.3) in [46.] • 



This can be identified as the Newtonian energy equa- 
tion for the coupled fluid. 

Equations (jA14[) and (jA17p can be expressed in more 
familiar forms by introducing the total (Lagrangian) time 
derivative comoving with the fluid: 



d 

d7 



di 



(A18) 



The continuity equation V^(pow^) = gives dpo/dr = 
—poV^u^. Hence, 



Po dr 



(A19) 



Combining Eqs. dZHI, C^TSj) and (UT9|) yields 

E 



-^{poe + E) - —{poe 
dr Pa 

+V^F^ + F^a^ = , 



dr 



(A20) 

Dividing both sides of Eq. (|A20p by po a-nd writing Ftot = 
Poe + E and Ftot = F + P, we get 



_d_ 

d^ 



E,. 



-Pt. 



1 



1 

PO 



V^F'^ 



1 



F^c 



Po / dr \po/ Po ' Po 

(A21) 

In Newtonian limit, d/dr d/dt, and we may neglect 
F'^Qfi and 9fF° as explained above, and Eq. (jA2ip re- 
duces to 



d_ /Ftot 
dt \ Po 



-Ft, 



1 

dt \po 



-V • F 



PO 



(A22) 



which is the familiar first-law of thermodynamics in the 
case where entropy is generated by radiation. 



d. Radiation equations 
From Eqs. §^ and we get 



dtf + d.,{a^^R 



iO\ 



G,j (A23) 
(A24) 



APPENDIX B: ESTIMATION OF 
CHARACTERISTIC SPEEDS 

In order to compute HLL fluxes, we must compute 
the maximum left-going wave speed c_ and maximum 
right-going wave speed c+ on both sides of the interface. 
We estimate the wave speed by computing the dispersion 
relation due to a small perturbation on a magnetized, ra- 
diating plasma of uniform po, P, , E and F*. In the 
comoving frame, = 0. We further choose our coordi- 
nate system so that the spacetime is locally Minkowski 
(i.e. g^i, = rj^^). To compute the dispersion relation, we 
consider a perturbation of the form 



Po 


= Po - 






p 


= FH 




t) 

7 


u' 


= Su' 


i(fecm-iC-Wcmt) 
^ •> 




B' 


= B' 






E 


= FH 


^ ,5Fe'(^"'-^-""' 


t) 

7 


pi 


= F' 







(Bl) 



where bar denotes the unperturbed quantity, and the 
subscript "cm" refers to comoving frame values. Substi- 
tuting these 12 expressions into our 12 radiation-MHD 
equations dSS]), jM]), gH), ^ and ^ and keeping 
terms linear in the perturbation, we obtain a matrix 
equation of the form MX = 0. Here iW is a 12 x 12 
matrix, and X = {5po SP Su' SB' 5E SF'f, 
where superscript t denotes the transpose. For simplic- 
ity, we drop the radiative source terms in deriving 
the dispersion relation (the curvature source terms van- 
ish in Minkowski spacetime). The dispersion relation is 
obtained by setting det(A4') — 0, which leads, after some 
algebra, to the following equation: 



4 / 

where 



2 

cm 



^cm/3)Km - (fccm ' V a) ]Q(t^cm, fccm) 



cmi ^cm) 



[k. 

's i^c 



2 

cm m 



VA?. 



= , 
(B2) 



(B3) 
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Here Cg = ^/TP/poh is the sound speed, va = 
\/h'^ I {poh + is the Alfven speed, and Cm = 
\/v\ + Cs(l — v\). The solution ujcm — corresponds to 
pure, stationary density perturbations, uj'^^ — is 
related to the propagation speed of the radiation flux for 
the nearly isotropic radiative diffusion, uj'^^ = (fccm •'I'a)^ 
corresponds to the Alfven waves, and QiiVcm, fccm) = 
corresponds to magnetosonic waves. As in j54{, we re- 
place the dispersion relation (5(i^cmi fecm) = by uj'^^ — 
Cm^cm = as it is more convenient when calculating the 
characteristic speed in the grid frame. As pointed out 
in [13], this modified dispersion relation overestimates 
the maximum wave speed by a factor of < 2 in the co- 
moving frame. 

Since the HLL scheme only requires the information 
on the maximum and minimum characteristic speeds, we 
use the following dispersion relation to estimate the char- 
acteristic speeds: 



(B4) 



APPENDIX C: ANALYTIC SOLUTIONS FOR 
RADIATION TESTS IN MINKOWSKI 
SPACETIME 

1. One-Dimensional waves in Minkowski Spacetime 




U2, and C/3. In the presence of a shock, across which 
P is discontinuous, these "conserved quantities" give the 
Rankine-Hugoniot conditions, relating P's on both sides 
of the shock. The remaining two ODEs [Eqs. (|C4p and 
(|C5p ] can be integrated numerically using, for example, 
a fourth-order Runge-Kutta integrator. During the inte- 
gration, we need to compute P from U, which we outline 
as follows. For simplicity, we consider the case without 
a large-scale electromagnetic field {T^J^^ = — B^). Us- 
ing Eq. dST]) for Tt"", Eq. ^ for i?'^'', and combining 
Eqs. dUIll-jCl]), we obtain 



[po + {n+l)P] u^'u" = U2-U4 
+ P = U.-U. 



[po + {n+l)P] {i 



= Ua 



3 



-E 



(C6) 
(C7) 
(C8) 

(C9) 



= C/4 

= U5 .(CIO) 



Eliminating po and P from Eqs. 
expression for u^: 



C6])-(IC8]), we get an 



Uiu" + {n+ l)Ubu"u^ - Ua[{u°f + n(u^)2] = (Cll) 



where u*^ = -^/l + (u^)^. Hence can be determined by 
solving the above algebraic equation. Having obtained 
u^, the quantities po and P are computed from 



We begin by assuming a F-law EOS and write T = 
1 + 1/n, where n is the polytropic index. Hence the EOS 
becomes poe = nP. 

We consider a stationary, infinite fluid in Minkowski 
spacetime with planar symmetry, hence we drop all time 
derivatives, all y and z components, and all y and z 
derivatives. It follows from Vp(pou^) = 0, S/^Tf"" = 
and \7^R^"' = that 





= , 


(CI) 


rpOx 


= 0, 


(C2) 


rpXX 


= 0, 


(C3) 




= -G" , 


(C4) 




= -G^ . 


(C5) 



It is convenient to define 



( Po\ 








( ° \ 


p 













,u = 




and S = 





E 




j^Ox 






\f- J 








v-G-y 



The system of ODEs in ([CT|) - ([C5|) becomes 

dxViP) = s(P) . 

The first three equations [Eqs. (|Cip - (|C3p ] are read- 
ily integrated, giving three "conserved quantities" Ui, 



Po 
P 



Ui 



Ua- 



(C12) 
(C13) 



The values of E and F^ are obtained by solving the linear 
Eqs. (IC9| and (|ClOl) . The result is 



E 

px 



A 
A 



where 



A, 



Af 



3 
4 



^ -fl 



U5 



(C14) 
(C15) 

(C16) 
(C17) 

(C18) 



We note that Eq. (jClip will, in general, have more 
than one real root. Indeed, this must be the case in 
order for shocks to exist. In the absence of a shock, the 
appropriate root is chosen by continuity. 

Our one-dimensional tests can be divided into two 
groups: fully continuous configurations, and discontin- 
uous configurations (i.e. shocks are present). In either 
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case, we begin by specifying boundary conditions on the 
asymptotic left side (a; = — oo). In practice, we set up 
our computational domain with x G [— i, L]. We specify 
POi Pi a-iid at X = —L, denoting them as poL, Pl and 
We also impose that the radiation be in LTE with 
the gas a,i X — —L [El = cirTI = aiim'^{PL/ pqlY] and 
set _F£ = /i, where is a parameter chosen so that 
Ff^/EL ^ 1. Thus, we have specified 5 boundary condi- 
tions for our 5 ODE's. The values of J7i, U2 and C/3, which 
are independent of are determined. If all quantities in 
the configuration are continuous throughout the compu- 
tational domain, these boundary conditions at a; = —L 
arc sufficient for us to integrate Eqs. (jC9[) and (|C10[) from 
X = —L to X = L. 

If however, we wish to determine a configuration that 
contains a discontinuity, which we set at x = 0, we use a 
shooting method, described as follows. In addition to the 
boundary conditions set at a: = —L, we also demand that 
the radiation be in LTE with the gas, and that F^ — fji 
ai X — L. The parameter fn is chosen so that F'-' /E <C 1 
at a; = L. The constants Ui, U2 and U3 are fixed, giving 
the values of P at a; = — L. Denote the values of po, 
P, and i? at a; = L by pofl, Pr, and w^, and Er, 
respectively. The LTE condition at a; = L gives Er = 
aRm'^{PR/ Por)'^. From the definition of [/i, U2 and t/s, 
it is straightforward to show that 



POR 



-fRU^R 



U2 - Uiul + 4 [(72 - C/34"fl] 

[{n ~ 3)C/i<] , (C19) 



then integrate Eqs. I|C4[) and (jC5|) from both x = ±L to 
X = 0. Since J7i, U2 and U-j, are constants, the Rankine- 
Hugoniot junction conditions are automatically satisfied 
at the shock front (a; — 0). Hence, we only have to im- 
pose the junction conditions for the radiation variables, 
which are the continuity of _R°^ and i?^^. In the New- 
tonian limit, these conditions reduce to the continuity of 
E and F^ at the shock front, but this is not the case in 
general. In any case, we need two junction conditions 
at the shock front, so a well-posed shooting problem re- 
quires varying two boundary condition parameters until 
the solution satisfies the two junction conditions at the 
shock front. We use /l and Jr as such two parameters. 

For a pure hydrodynamic (or MHD) shock, the profiles 
of P are constants on each side of the shock front. This is 
not the case for a radiating hydrodynamic shock, where 
P vary with x and approach constants only in the asymp- 
totic regions {x — s- ±00). This variation results from the 
radiative source terms and , given by Eq. ([33]), 
which vanishes when the radiation and fluid are in strict 
LTE and the radiation flux vanishes (i.e. in the asymp- 
totic regions). The length scale over which the parame- 
ters vary between the two asymptotic regions is a few op- 
tical depths. We need to choose k to ensure that x = ±L 
are in the asymptotic regions. In practice, one can esti- 
mate the optical depth r 
and choose k so that r » 1 
T ^ 10. Choosing a larger k does not change the proflle 
(plotted against the rescaled coordinate kx) significantly. 



2. Special Analytic Case 



'j^poKdx ~ (poL -\- Por.)kL 
For our tests, we choose 



Uiu% 
+ fR 



Pr 

POR 
2 





fPR\ 


+ ^aRin 






KporJ 



Ur>U 



R"'R 



U2 = 



(C20) 



where 



I = y/TTJv^. Substituting Eq. ((CT9l) into 

Eq. (|C20[) gives an algebraic equation for u'^, which can 
be solved numerically. Obviously, wfj. = ii| is a solution, 
but we look for another solution in order to obtain a 
configuration containing a shock. Having determined u^, 
the other quantities are computed as follows: 



POR 

Pr 



Ul/U^R 

Pr 



POR 



-POR 



Er = aRiJi 



{- 

\POR 



(C21) 
(C22) 

(C23) 



where Pr/pqr in above equations are computed from 
Eq. (jCigp . To generate the ID shock configuration, we 
first specify p^, P^, u^, and consider and Jr as free 
parameters. The other quantities at a; = ±L are fixed by 
the LTE condition at a; = ±L and Eqs. ((CT9)) - ((C23) ). We 



In Newtonian limit, the solutions of (|C6|) - (|C10[) can 
be written in analytic form under special conditions, as 
stated in [46, 64]. Here we briefly summarize this solu- 
tion. 

In Newtonian limit, Eqs. (|C6|) -- (|C10p become 



Pov 



Ui (C24) 

poe + P + E + V]v + F'' = U2 (C25) 



pov^ + P + V = C/3 ,(C26) 

while by dropping time derivatives in Eqs. (jA29|l - (|A30[l . 
we get 



px 



AnB - E 

1 dE(r) 

~3 dr 



(C27) 
(C28) 



where v = v'^ and r is the optical depth, given by 
dr = PQKdx. As in [1^, HH, we consider a strong shock 
propagating into cold gas, so that the pressure and inter- 
nal energy of the unshocked gas (a: < 0) can be neglected. 
We also assume that the gas is optically thick so that we 
may use the diffusion approximation, and that it is suf- 
ficient to account for the radiation fiux, while neglecting 
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radiation energy density and radiation pressure. Under 
these assumptions, the above equations can be rewritten 
as: 



Pov ^ Pqlvl 
pov^ + P = PQLvl 

,2 /n\ , IPX _ „ „ ,2 



Pov{poe + v'/2) + F''' = poLvi/2 
Combining Eqs. ((U28)) and ([UST]) gives 



(C29) 
(C30) 
(C31) 



(C32) 



dr^ ' ^ dr 

Solvin g th ese coupled equations using methods outlined 
in [46l. |64| one arrives at: 
Preshock Medium {x < 0): 

1 



2V3 



1 



Postshock Medium {x > 0): 



-V3|r| 



La^T^e-^l^l 

E ^ aflT4(l - l/2e-^l^l) 



(C33) 
(C34) 

(C35) 
(C36) 



Here Tr is the asymptotic temperature in the postshock 
region and r is measured from the shock front (i.e. t{x) — 
Kp{x')dx'). 



APPENDIX D: THERMAL 
OPPENHEIMER-SNYDER SOLUTION 

Here we summarize the analytic solutions derived 
in [2^ for a strong field (black hole) dynamical sce- 
nario, used to compare with our numerical results in Sec- 
tion [iVBl 

The standard Oppenheimer-Snyder (OS) collapse so- 
lution for a homogeneous dust ball was first derived 
In [2^, the collapsing sphere is subjected to 
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radiation and gas pressure perturbations which are as- 
sumed to be sufficiently small that the spacetime metric 
and density evolution are well-approximated by the OS 
solution. For this to be true, we require P/po M/R 
and V / po = E/3po ^ M/R for a star with mass M 
and radius R. We are interested in the evolution of the 
radiation quantities E and F^ inside the star. 

Inside the star, the OS metric is given by the closed 
Friedmann line element 

ds^ = -dr^ + a^{T){dx^ + sin^ xd^^) (Dl) 

where r is the proper time of a fluid element and x is 
a Lagrangian radial coordinate. The scale factor a(T) is 
given in parametric form according to 

a = ^am(l + cosry) , (D2) 



a™(?? -t- sin 77) 



(D3) 



Here 77 is the conformal time and constant which 

is related to the initial areal radius Ri of the star. (The 
subscript i denotes initial values.) The radius is given by 



R = ^i?i(l + cos?7) . 



The exterior Schwarzschild line element is 

-1 



(D4) 



ds^ 



1 ] dt^ 



(D5) 

where Ts is the Schwarzschild (areal) radius. Matching 
the two metrics at = i? gives 



2M 



(D6) 



Denote xo ^ the Lagrangian radial coordinate x the 
stellar surface. It follows from Eqs. (|DT|1 . ((D5)) . 02\i and 
(lD4l) that 



R 2M 



a y R, 



(D7) 



In Friedmann comoving coordinates the density po is al- 
ways homogeneous and given by 



Pq_ 

POi 



Q- 



where 



a 1 
Q = — = x(l + cosr/) 

Clrn ^ 



(D8) 



(D9) 



An analytic solution for the E and F is derived in [29|, 
Isoj by assuming (1) diffusion approximation, (2) that the 
radiation and fluid are in LTE {E = ajiT"^), and (3) 
that radiation pressure is much greater than gas pressure 
{V ^ P). We summarize the solution below. 

Define the radiation energy density and flux "cor- 
rected" for adiabatic contraction as Ec = Q*E and 
Fc — Q'^F respectively. Define a time parameter f as 



1 



Ri ( sinxo 



4 . 1 . „ 

, , , , I 1 7? H — sm 77 H — sm 277 

4(T« + r«)V8MV Xo y V 3 ' 6 ' 

(DIO) 

where r" and are the initial absorption and scattering 
optical depths respectively, related to the absorption and 
scattering opacities and by r° = K^poi-Ri and = 
K'^poiRi. (Note that the analytic solution assumes k° 
and to be constant throughout the collapse.) Define 
also a normalized Lagrangian radius z = x/xoi such that 
< z < 1 within the star. Then the interior corrected 
energy density is given by 



E,if,z) = 2E, 



sin{xQz] 



0X0 



sin(n7rz) 



77 TT 



(-1) 



n+1 



7727r2 - xl 



(Dll) 
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where Ei is the initial value of E, assumed to be con- 
stant throughout the star. Note that in the derivation 
of Eq. (jPTTIl the "zero temperature" boundary condition 
{E = a.t the stellar surface) has been used. 

The diffusion approximation gives the expression for 
the corrected radiation flux Fc in terms of the gradient 
of 



Fc{f,z) 



2 E.Q'^ 1 / sinxo 



(D12) 



3 T" + Xo Vsin(xoz) 



e " '^[xo sin(7i7rz) cos(xo-2) 
— UTT sin(xo2) cos(n7rz)] 



(D13) 



Finally, an analytic expression for the ideal gas pres- 
sure P — pqT /m is obtained from the LTE condition 



E = grT^ 



aRin 



Hence 



E 



1/4 



(DM) 



where po and E are given by the analytic solutions above. 

A comparison between the analytic solution for ther- 
mal OS collapse and the exact solution of the Boltzmann 
equation of radiative transfer for the same problem is 
given in (30j . 



[1] M. C. Begelman, R. D. Blandford, and M. J. Rees, Re- 
views of Modern Physics 56, 255 (1984). 

[2] T. Piran, Reviews of Modern Physics 76, 1143 (2005), 
arXiv:astro-ph/0405503. 

[3] S. A. Balbus and J. F. Hawley, Reviews of Modern 
Physics 70, 1 (1998). 

[4] J. M. Miller, J. Raymond, A. Fabian, D. Steeghs, 
J. Homan, C. Reynolds, M. van der Klis, and R. Wij- 
nands. Nature (London) 441, 953 (2006), arXiv:astro- 
ph/0605390. 

[5] T. W. Baumgarte, S. L. Shapiro, and M. Shibata, Astro- 

phys. J. Lett. 528, L29 (2000), arXiv:astro-ph/9910565. 
[6] S. L. Shapiro, Astrophys. J. 544, 397 (2000), 

arXiv:astro-ph/0010493. 
[7] M. D. Duez, Y. T. Liu, S. L. Shapiro, M. Shibata, and 

B. C. Stephens, Physical Review Letters 96, 031101 

(2006), arXiv:astro-ph/05 10653. 
[8] M. D. Duez, Y. T. Liu, S. L. Shapiro, M. Shibata, 

and B. C. Stephens, Phys. Rev. D 73, 104015 (2006), 

arXiv:astro-ph/0605331. 
[9] H.-T. Janka, K. Langanke, A. Marek, G. Martmez- 

Pinedo, and B. Midler, Phys. Rep. 442, 38 (2007), 

arXiv ;astro-ph /06 1 2072 . 
[10] A. Burrows, L. Dessart, C. D. Ott, and E. Livne, Phys. 

Rep. 442, 23 (2007), arXiv;astro-ph/0612460. 
[11] A. Mezzacappa, S. W. Bruenn, J. M. Blondin, W. R. 

Hix, and O. E. Bronson Messer, in American Institute of 

Physics Conference Series (2007), vol. 924 of American 

Institute of Physics Conference Series, pp. 234-242. 
[12] D. Richstone, E. A. Ajhar, R. Bender, G. Bower, 

A. Dressier, S. M. Faber, A. V. Filippenko, K. Gebhardt, 

R. Green, L. C. Ho, et al.. Nature (London) 395, A14-I- 

(1998), arXiv:astro-ph/9810378. 
[13] L. Ho, in Observational Evidence for the Black Holes in 

the Universe, edited by S. K. Chakrabarti (1999), vol. 234 

of Astrophysics and Space Science Library, pp. 157 — h. 
[14] T. W. Baumgarte and S. L. Shapiro, Astrophys. J. 526, 

941 (1999), arXiv:astro-ph/9909237. 
[15] S. L. Shapiro, in Coevolution of Black Holes and Calaxies, 

edited by L. C. Ho (2004), pp. 103— 
[16] S. L. Shapiro, Astrophys. J. 620, 59 (2005), arXiviastro- 

ph/0411156. 



[17] M. D. Duez, Y. T. Liu, S. L. Shapiro, and B. C. 
Stephens, Phys. Rev. D 72, 024028 (2005), arXiv:astro- 
ph/0503420. 

[18] M. Shibata and T. Nakamura, Phys. Rev. D 52, 5428 
(1995). 

[19] T. W. Baumgarte and S. L. Shapiro, Physical Review D 

59, 024007 (1999). 
[20] G. Toth, Journal of Computational Physics 161, 605 

(2000). 

[21] M. Shibata, Y. T. Liu, S. L. Shapiro, and B. C. 
Stephens, Phys. Rev. D 74, 104026 (2006), arXiv:astro- 
ph/0610840. 

[22] Y. T. Liu, S. L. Shapiro, and B. C. Stephens, Phys. Rev. 

D 76, 084017 (2007), arXiv:0706.2360. 
[23] Z. B. Etienne, J. A. Faber, Y. T. Liu, S. L. Shapiro, 

and T. W. Baumgarte, Phys. Rev. D 76, 101503 (2007), 

arXiv:0707.2083. 
[24] Z.B. Etienne, J.A. Faber, Y.T. Liu, S.L. Shapiro, 

K. Taniguchi, and T.W. Baumgarte, submitted to Phys. 

Rev. D (arXiv:0712.2460). 
[25] R. W. Lindquist, Annals of Physics 37, 487 (1966). 
[26] K. S. Thorne, Mon. Not. R. Astron. Soc. 194, 439 (1981). 
[27] J. Schmid-Burgk, Astrophysics and Space Science 56, 

191 (1978). 

[28] P. J. Schinder, Phys. Rev. D 38, 1673 (1988). 
[29] S. L. Shapiro, Phys. Rev. D 40, 1858 (1989). 
[30] S. L. Shapiro, Astrophys. J. 472, 308 (1996). 
[31] P. J. Schinder and S. A. Bludman, Astrophys. J. 346, 
350 (1989). 

[32] A. Mezzacappa and R. A. Matzner, Astrophys. J. 343, 
853 (1989). 

[33] S. Zane, R. Turolla, L. Nobih, and M. Erna, Astrophys. 
J. 466, 871 (1996), arXiv:astro-ph/9602034. 

[34] M. Liebendorfer, O. E. B. Messer, A. Mezzacappa, S. W. 
Bruenn, C. Y. Cardall, and F.-K. Thielemann, Astro- 
phys. J. Supp. 150, 263 (2004), arXiv:astro-ph/0207036. 

[35] K. S. Thorne, R. A. Flammang, and A. N. Zytkow, Mon. 
Not. R. Astron. Soc. 194, 475 (1981). 

[36] L. RezzoUa and J. C. Miller, Classical and Quantum 
Gravity 11, 1815 (1994), arXiv:astro-ph/9406055. 

[37] L. Zampieri, ,L C. Miller, and R. Turolla, Mon. Not. R. 
Astron. Soc. 281, 1183 (1996), arXiv:astro-ph/9607030. 



20 



[38] R. TuroUa, S. Zano, L. Zampieri, and L. Nobili, Mon. 
Not. R. Astron. Soc. 283, 881 (1996), arXiv:astro- 
ph/9608174. 

[39] L. Rezzolla and J. C. Miller, Phys. Rev. D 53, 5411 

(1996), arXiv:astro-ph/9510039. 
[40] S. Balborg, L. Zampieri, and S. L. Shapiro, Astrophys. J. 

541, 860 (2000), arXiv:astro-ph/0004234. 
[41] M. D. Ducz, R Marronetti, S. L. Shapiro, and T. W. 

Baurngarte, Physical Review D 67, 024004 (2003). 
[42] M. Alcubierre, B. Briigmann, D. Polhicy, E. Scidel, and 

R. Takahashi, Phys. Rev. D 64, 061501 (2001), arXiv:gr- 

qc/0104020. 

[43] M. D. Duez, S. L. Shapiro, and H.-J. Yo, Physical Review 

D 69, 104016 (2004). 
[44] M. CampaneUi, C. O. Lousto, P. Marronetti, and Y. Zlo- 

chower. Physical Review Letters 96, 111101 (2006), 

arXiv:gr-qc/0511048. 
[45] J. G. Baker, J. CentreUa, D.-I. Choi, M. Koppitz, and 

J. van Meter, Physical Review Letters 96, 111102 (2006), 

arXiv:gr-qc/0511103. 
[46] D. Mihalas and B. Weibel Mihalas, Foundations of radia- 
tion hydrvdynamics (New York: Oxford University Press, 

1984, 1984). 

[47] M. Alcubierre, B. Briigmann, D. Holz, R. Takahashi, 
S. Brandt, E. Seidel, J. Thornburg, and A. Ashtekar, In- 
ternational Journal of Modern Physics D 10, 273 (2001), 
arXiv:gr-qc/9908012. 

[48] M. D. Ducz, Y. T. Liu, S. L. Shapiro, and B. C. 
Stephens, Phys. Rev. D 69, 104030 (2004), arXiv:astro- 
ph/0402502. 

[49] http://www.cactuscode.org/. 

[50] M. D. Ducz, S. L. Shapiro, and H.-J. Yo, Physical Review 

D 69, 104016 (2004). 
[51] B. van Leer, Journal of Computational Physics 23, 276 

(1977). 

[52] A. Marten, P. D. Lax, and v. B. J., SIAM Rev. 25, 35 
(1983). 

[53] C. R. Evans and J. F. Hawley, Astrophys. J. 332, 659 



(1988). 

[54] C. F. Gammic, J. C. McKinncy, and G. Toth, Astrophys. 

J. 589, 444 (2003), arXiv:astro-ph/0301509. 
[55] T. W. Baumgarte and S. L. Shapiro, Astrophys. J. 585, 

930 (2003), arXiv:astro-ph/0211339. 
[56] M. Rampp and H.-T. Janka, Astronomy and Astro- 
physics 396, 361 (2002), arXiv:astro-ph/0203101. 
[57] M. Liebendorfer, A. Mezzacappa, F.-K. Thielemann, 

O. E. Messer, W. R. Mix, and S. W. Bruenn, Phys. Rev. 

D 63, 103004 (2001), arXiv:astro-ph/0006418. 
[58] A. Burrows, T. Young, P. Pinto, R. Eastman, and T. A. 

Thompson, Astrophys. J. 539, 865 (2000), arXiv:astro- 

ph/9905132. 

[59] C. D. Levermore and G. C. Pomraning, Astrophys. J. 

248, 321 (1981). 
[60] M. Liebendoerfer, S. C. Whitehouse, and T. Fischer, 

ArXiv e-prints 711 (2007), 0711.2929. 
[61] J. A. Font, M. Miller, W.-M. Suen, and M. Tobias, Phys. 

Rev. D 61, 044011 (2000), arXiv:gr-qc/9811015. 
[62] M. Shibata and Y.-I. Sckiguchi, Phys. Rev. D 72, 044014 

(2005), arXiv:astro-ph/0507383. 
[63] J. A. Faber, T. W. Baumgarte, Z. B. Eticrmc, S. L. 

Shapiro, and K. Taniguchi, Phys. Rev. D 76, 104021 

(2007), arXiv:0708.2436. 
[64] Y. B. Zel'Dovich and Y. P. Raizer, Physics of shock waves 

and high-temperature hydrodynamic phenomena (New 

York: Academic Press, 1966/1967, edited by Hayes, 

W.D.; Probstein, Ronald F., 1967). 
[65] J. R. Oppenheimer and H. Snyder, Physical Review 56, 

455 (1939). 

[66] L. I. Petrich, S. L. Shapiro, and S. A. Teukolsky, Phys. 
Rev. D 31, 2459 (1985). 

[67] M. Shibata, T. W. Baumgarte, and S. L. Shapiro, Phys- 
ical Review D 61, 044012 (2000). 

[68] I. D. Novikov and K. S. Thorne, in Black Holes, edited 
by C. DeWitt and B. S. DeWitt (Gordon and Breach, 
New York, 1973), p. 343. 



